Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(brms)
library(broom.mixed)
library(tidybayes)
library(bayesplot)
library(standist) # for visualizing distributions
library(rstanarm)
library(ggeffects)
library(rstan)
library(DHARMa)
library(ggridges)
source("helperFunctions.R")

Scenario

Crab_shrimp_coral

To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a predatory sea star and one of four symbiont combinations:

no symbiont,
a crab symbiont
a shrimp symbiont
both a crab and shrimp symbiont.

The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, sea stars and symbiont.

The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.

Read in the data

mckeon <- read_csv("../public/data/mckeon.csv", trim_ws = TRUE)
mckeon %>% glimpse()

## Rows: 80
## Columns: 3
## $ BLOCK     <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10…
## $ PREDATION <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ SYMBIONT  <chr> "none", "none", "none", "none", "none", "none", "none", "non…

Since the response here is the presence or absence of predation (feeding scars), a binomial distribution is appropriate.

We need to make sure that the categorical variable and the random effect are defined as factors. When doing so, it might be valuable to rearrange the order of the fixed effect (SYMBIONT) such that the none group is considered the first group. This way, the other levels will all naturally be compared to this level (hence it will be treated as a reference of control group).

mckeon <- mckeon %>%
  mutate(
    BLOCK = factor(BLOCK),
    SYMBIONT = factor(SYMBIONT, levels = c("none", "crabs", "shrimp", "both"))
  )

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are vectors of the fixed and random effects parameters respectively and $\bf{X}$ is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation. $\bf{Z}$ represents a cell means model matrix for the random intercepts associated with individual coral colonies.

ggplot(mckeon, aes(y = PREDATION, x = SYMBIONT)) +
  geom_point(position = position_jitter(width = 0.2, height = 0)) +
  facet_wrap(~BLOCK)

Fit the model

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

mckeon.rstanarm <- stan_glmer(PREDATION ~ SYMBIONT + (1 | BLOCK),
  data = mckeon,
  family = binomial(link = "logit"),
  iter = 5000,
  warmup = 2000,
  chains = 3,
  thin = 5,
  refresh = 0,
  cores = 3
)

mckeon.rstanarm %>% prior_summary()

## Priors for model '.' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
##   Adjusted prior:
##     ~ normal(location = [0,0,0], scale = [5.74,5.74,5.74])
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

This tells us:

for the intercept (when the family is binomial), a normal prior with a mean of 0 and a standard deviation of 2.5 is used. These are the default priors for bernoulli models. A mean of 0 on a logit scale corresponds to a probability of 0.5. Hence the prior expected value is 0.5.
for the coefficients (in this case, just the slope), the default prior is a normal prior centred around 0 with a standard deviation of 2.5. For binomial models, this is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the predictor (then rounded).

model.matrix(~SYMBIONT, data = mckeon) %>%
  apply(2, sd) %>%
  (function(x) 2.5 / x)

##    (Intercept)  SYMBIONTcrabs SYMBIONTshrimp   SYMBIONTboth 
##            Inf       5.737305       5.737305       5.737305

Assessing priors

Lets now run with priors only so that we can explore the range of values they allow in the posteriors.

mckeon.rstanarm1 <- update(mckeon.rstanarm, prior_PD = TRUE)

mckeon.rstanarm1 %>%
  ggpredict() %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

## $SYMBIONT

Conclusions:

it is very difficult to assess the priors from the predicted posteriors when the posteriors are on the response scale as they will always approach values of 0 and 1 (no matter how wide they are on the link scale).

Interestingly, if we instead use ggemmeans, it will (perhaps erroneously) generate the partial effects on the link scale (yet with an incorrectly labelled y-axis). Probability scale values of 0.99 and 0.01 (very large and small respectively) correspond to value of approximately 4.6 and -4.6 respectively on the logit scale.

In the following partial plot, the routine attempts to convert the y-axes tick marks into percentages. Hence a value of 0.1 is converted to 10%. Consequently, values of -5 and 5 are labelled as -500% and 500% respectively. We know that on the probability scale, values must asymptote towards 0 and 1. Posterior intervals beyond -5 and 5 might be considered wide.

mckeon.rstanarm1 %>%
  ggemmeans(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

Alternatively, we could create the plot ourselves…

mckeon.rstanarm1 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  as.data.frame() %>%
  ggplot(aes(y = emmean, x = SYMBIONT)) +
  geom_hline(yintercept = c(5, -5), linetype = "dashed") +
  geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
  geom_point(
    data = mckeon, aes(y = PREDATION),
    position = position_jitter(width = 0.2, height = 0),
    alpha = 0.4, color = "red"
  )

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

$\beta_0$: normal centred at 0 (since this is 0.5 on the logit scale) with a standard deviation of 2
- mean of 0: since logit(0.5)
- alternatively, an argument could be made for a mean of 0.51: since logit(mean(mckeon$PREDATION))
- sd of 2: since this not too large on logit scale
$\beta_1$: normal centred at 0 with a standard deviation of 0.5
- sd of 0.5: since model.matrix(~SYMBIONT, data=mckeon)[,-1] %>% apply(2,sd)
$\Sigma$: decov with:
- regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

I will also overlay the raw data for comparison.

mckeon.rstanarm2 <- stan_glmer(PREDATION ~ SYMBIONT + (1 | BLOCK),
  data = mckeon,
  family = binomial(link = "logit"),
  prior_intercept = normal(0, 2, autoscale = FALSE),
  prior = normal(0, 1, autoscale = FALSE),
  prior_covariance = decov(1, 1, 1, 1),
  prior_PD = TRUE,
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

mckeon.rstanarm2 %>%
  ggpredict(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

mckeon.rstanarm2 %>%
  ggemmeans(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

Now lets refit, conditioning on the data.

mckeon.rstanarm3 <- update(mckeon.rstanarm2, prior_PD = FALSE)

Plotting prior and posterior

posterior_vs_prior(mckeon.rstanarm3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Conclusions:

in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggemmeans(mckeon.rstanarm3, ~SYMBIONT) %>% plot(add.data = TRUE)

ggpredict(mckeon.rstanarm3, ~SYMBIONT) %>% plot(add.data = TRUE)

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

mckeon.form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1|BLOCK),
                  family=binomial(link='logit'))
options(width=150)
mckeon.form %>% get_prior(data = mckeon)

##                 prior     class           coef group resp dpar nlpar bound       source
##                (flat)         b                                                 default
##                (flat)         b   SYMBIONTboth                             (vectorized)
##                (flat)         b  SYMBIONTcrabs                             (vectorized)
##                (flat)         b SYMBIONTshrimp                             (vectorized)
##  student_t(3, 0, 2.5) Intercept                                                 default
##  student_t(3, 0, 2.5)        sd                                                 default
##  student_t(3, 0, 2.5)        sd                BLOCK                       (vectorized)
##  student_t(3, 0, 2.5)        sd      Intercept BLOCK                       (vectorized)

options(width=80)

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

$\beta_0$: normal centred at 0 (since this is 0.5 on the logit scale) with a standard deviation of 2
- mean of 0: since logit(0.5)
- alternatively, an argument could be made for a mean of 0.51: since logit(mean(mckeon$PREDATION))
- sd of 2: since this not too large on logit scale
$\beta_1$: normal centred at 0 with a standard deviation of 0.5
- sd of 0.5: since model.matrix(~SYMBIONT, data=mckeon)[,-1] %>% apply(2,sd)
$\sigma_j$: half-cauchy with parameters 0 and 1.
$\Sigma$: decov with:
- regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large $sigma$ in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

half-t: as the degrees of freedom approach infinity, this will approach a half-normal
half-cauchy: this is essentially a half-t with 1 degree of freedom
exponential

I will also overlay the raw data for comparison.

mckeon %>%
  group_by(SYMBIONT) %>%
  summarise(
    mean(PREDATION),
    var(PREDATION)
  )

standist::visualize("normal(0,2)", xlim = c(0, 200))

standist::visualize("student_t(3, 0, 2.5)",
  "gamma(2,0.5)",
  "cauchy(0,1)",
  xlim = c(-10, 25)
)

priors <- prior(normal(0, 2), class = "Intercept") +
  prior(normal(0, 2), class = "b") +
  prior(cauchy(0, 1), class = "sd")

mckeon.form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1 | BLOCK),
  family = binomial(link = "logit")
)
mckeon.brm2 <- brm(mckeon.form,
  data = mckeon,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0
)

mckeon.brm2 %>%
  ggpredict(~SYMBIONT) %>%
  plot(add.data = TRUE)

mckeon.brm2 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  as.data.frame() %>%
  ggplot(aes(y = emmean, x = SYMBIONT)) +
  geom_hline(yintercept = c(5, -5), linetype = "dashed") +
  geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
  geom_point(
    data = mckeon, aes(y = PREDATION),
    position = position_jitter(width = 0.2, height = 0),
    alpha = 0.4, color = "red"
  )

Now lets refit, conditioning on the data.

mckeon.brm3 <- update(mckeon.brm2,
  sample_prior = "yes",
  refresh = 0
)
save(mckeon.brm3, file = "../ws/testing/mckeon.brm3")

mckeon.brm3 %>%
  ggpredict(~SYMBIONT) %>%
  plot(add.data = TRUE)

mckeon.brm3 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  as.data.frame() %>%
  ggplot(aes(y = emmean, x = SYMBIONT)) +
  geom_hline(yintercept = c(5, -5), linetype = "dashed") +
  geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
  geom_point(
    data = mckeon, aes(y = PREDATION),
    position = position_jitter(width = 0.2, height = 0),
    alpha = 0.4, color = "red"
  )

Plotting prior and posterior

mckeon.brm3 %>% get_variables()

##  [1] "b_Intercept"           "b_SYMBIONTcrabs"       "b_SYMBIONTshrimp"     
##  [4] "b_SYMBIONTboth"        "sd_BLOCK__Intercept"   "r_BLOCK[1,Intercept]" 
##  [7] "r_BLOCK[2,Intercept]"  "r_BLOCK[3,Intercept]"  "r_BLOCK[4,Intercept]" 
## [10] "r_BLOCK[5,Intercept]"  "r_BLOCK[6,Intercept]"  "r_BLOCK[7,Intercept]" 
## [13] "r_BLOCK[8,Intercept]"  "r_BLOCK[9,Intercept]"  "r_BLOCK[10,Intercept]"
## [16] "prior_Intercept"       "prior_b"               "prior_sd_BLOCK"       
## [19] "lp__"                  "accept_stat__"         "stepsize__"           
## [22] "treedepth__"           "n_leapfrog__"          "divergent__"          
## [25] "energy__"

mckeon.brm3 %>%
  hypothesis("SYMBIONTcrabs=0") %>%
  plot()

mckeon.brm3 %>%
  hypothesis("SYMBIONTshrimp=0") %>%
  plot()

mckeon.brm3 %>% SUYR_prior_and_posterior()

Random intercept/slope model

While we are here, we might like to explore a random intercept/slope model

priors <- prior(normal(0, 2), class = "Intercept") +
  prior(normal(0, 2), class = "b") +
  prior(cauchy(0, 1), class = "sd") +
  prior(lkj_corr_cholesky(1), class = "L")

mckeon.form <- bf(PREDATION | trials(1) ~ SYMBIONT + (SYMBIONT | BLOCK),
  family = binomial(link = "logit")
)
mckeon.brm4 <- brm(mckeon.form,
  data = mckeon,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 1000,
  chains = 3,
  thin = 5,
  refresh = 0,
  control = list(adapt_delta = 0.99)
)

save(mckeon.brm4, file = "../ws/testing/mckeon.brm4")

mckeon.brm4 %>% get_variables()

##  [1] "b_Intercept"                             
##  [2] "b_SYMBIONTcrabs"                         
##  [3] "b_SYMBIONTshrimp"                        
##  [4] "b_SYMBIONTboth"                          
##  [5] "sd_BLOCK__Intercept"                     
##  [6] "sd_BLOCK__SYMBIONTcrabs"                 
##  [7] "sd_BLOCK__SYMBIONTshrimp"                
##  [8] "sd_BLOCK__SYMBIONTboth"                  
##  [9] "cor_BLOCK__Intercept__SYMBIONTcrabs"     
## [10] "cor_BLOCK__Intercept__SYMBIONTshrimp"    
## [11] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp"
## [12] "cor_BLOCK__Intercept__SYMBIONTboth"      
## [13] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth"  
## [14] "cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth" 
## [15] "r_BLOCK[1,Intercept]"                    
## [16] "r_BLOCK[2,Intercept]"                    
## [17] "r_BLOCK[3,Intercept]"                    
## [18] "r_BLOCK[4,Intercept]"                    
## [19] "r_BLOCK[5,Intercept]"                    
## [20] "r_BLOCK[6,Intercept]"                    
## [21] "r_BLOCK[7,Intercept]"                    
## [22] "r_BLOCK[8,Intercept]"                    
## [23] "r_BLOCK[9,Intercept]"                    
## [24] "r_BLOCK[10,Intercept]"                   
## [25] "r_BLOCK[1,SYMBIONTcrabs]"                
## [26] "r_BLOCK[2,SYMBIONTcrabs]"                
## [27] "r_BLOCK[3,SYMBIONTcrabs]"                
## [28] "r_BLOCK[4,SYMBIONTcrabs]"                
## [29] "r_BLOCK[5,SYMBIONTcrabs]"                
## [30] "r_BLOCK[6,SYMBIONTcrabs]"                
## [31] "r_BLOCK[7,SYMBIONTcrabs]"                
## [32] "r_BLOCK[8,SYMBIONTcrabs]"                
## [33] "r_BLOCK[9,SYMBIONTcrabs]"                
## [34] "r_BLOCK[10,SYMBIONTcrabs]"               
## [35] "r_BLOCK[1,SYMBIONTshrimp]"               
## [36] "r_BLOCK[2,SYMBIONTshrimp]"               
## [37] "r_BLOCK[3,SYMBIONTshrimp]"               
## [38] "r_BLOCK[4,SYMBIONTshrimp]"               
## [39] "r_BLOCK[5,SYMBIONTshrimp]"               
## [40] "r_BLOCK[6,SYMBIONTshrimp]"               
## [41] "r_BLOCK[7,SYMBIONTshrimp]"               
## [42] "r_BLOCK[8,SYMBIONTshrimp]"               
## [43] "r_BLOCK[9,SYMBIONTshrimp]"               
## [44] "r_BLOCK[10,SYMBIONTshrimp]"              
## [45] "r_BLOCK[1,SYMBIONTboth]"                 
## [46] "r_BLOCK[2,SYMBIONTboth]"                 
## [47] "r_BLOCK[3,SYMBIONTboth]"                 
## [48] "r_BLOCK[4,SYMBIONTboth]"                 
## [49] "r_BLOCK[5,SYMBIONTboth]"                 
## [50] "r_BLOCK[6,SYMBIONTboth]"                 
## [51] "r_BLOCK[7,SYMBIONTboth]"                 
## [52] "r_BLOCK[8,SYMBIONTboth]"                 
## [53] "r_BLOCK[9,SYMBIONTboth]"                 
## [54] "r_BLOCK[10,SYMBIONTboth]"                
## [55] "prior_Intercept"                         
## [56] "prior_b"                                 
## [57] "prior_sd_BLOCK"                          
## [58] "prior_cor_BLOCK"                         
## [59] "lp__"                                    
## [60] "accept_stat__"                           
## [61] "stepsize__"                              
## [62] "treedepth__"                             
## [63] "n_leapfrog__"                            
## [64] "divergent__"                             
## [65] "energy__"

mckeon.brm4 %>%
  hypothesis("SYMBIONTcrabs=0") %>%
  plot()

mckeon.brm4 %>%
  hypothesis("SYMBIONTshrimp=0") %>%
  plot()

mckeon.brm4 %>% SUYR_prior_and_posterior()

mckeon.brm4 %>%
  posterior_samples() %>%
  dplyr::select(-`lp__`) %>%
  pivot_longer(everything(), names_to = "key") %>%
  filter(!str_detect(key, "^r")) %>%
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    ## Class = ifelse(str_detect(key, 'Intercept'),  'Intercept',
    ##         ifelse(str_detect(key, 'b'),  'b', 'sigma')),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_SYMBIONT.*|prior_b_SYMBIONT.*") &
        !str_detect(key, ".*:.*") ~ "SYMBIONT",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) %>%
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge()) +
  facet_wrap(~Class, scales = "free")

Compare models

We can compare the two models using LOO

(l.1 <- mckeon.brm3 %>% loo())

## 
## Computed from 2400 by 80 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo    -29.0  6.2
## p_loo         8.2  2.2
## looic        58.0 12.3
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     76    95.0%   1343      
##  (0.5, 0.7]   (ok)        1     1.2%   971       
##    (0.7, 1]   (bad)       3     3.8%   78        
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.

(l.2 <- mckeon.brm4 %>% loo())

## 
## Computed from 2400 by 80 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo    -23.2  5.1
## p_loo         9.8  2.6
## looic        46.5 10.2
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     59    73.8%   644       
##  (0.5, 0.7]   (ok)       17    21.2%   253       
##    (0.7, 1]   (bad)       4     5.0%   30        
##    (1, Inf)   (very bad)  0     0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.

loo_compare(l.1, l.2)

##   elpd_diff se_diff
## .  0.0       0.0   
## . -5.8       2.9

It appears that the random intercept/slope model is no better than the simpler random intercept model.

mckeon.brm4 %>%
  posterior_samples() %>%
  dplyr::select(-`lp__`) %>%
  pivot_longer(everything(), names_to = "key") %>%
  filter(!str_detect(key, "^r")) %>%
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_SYMBIONT.*|prior_b") ~ "TREATMENT",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) %>%
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
  facet_wrap(~Class, scales = "free")

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

pars <- mckeon.brm3 %>% get_variables()
pars <- pars %>%
  str_extract("^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*") %>%
  na.omit()
pars

## [1] "b_Intercept"         "b_SYMBIONTcrabs"     "b_SYMBIONTshrimp"   
## [4] "b_SYMBIONTboth"      "sd_BLOCK__Intercept"
## attr(,"na.action")
##  [1]  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## attr(,"class")
## [1] "omit"

mckeon.brm3 %>% mcmc_plot(type = "trace", variables = pars)

# OR
mckeon.brm3 %>% mcmc_plot(
  type = "trace",
  variable = "^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*",
  regex = TRUE
)

The chains appear well mixed and very similar

acf_bar (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

mckeon.brm3 %>% mcmc_plot(type = "acf_bar", variable = pars)

## OR
mckeon.brm3 %>% mcmc_plot(
  type = "acf_bar",
  variable = "^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*",
  regex = TRUE
)

There is no evidence of auto-correlation in the MCMC samples

rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mckeon.brm3 %>% mcmc_plot(type = "rhat_hist")

All Rhat values are below 1.05, suggesting the chains have converged.

neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

mckeon.brm3 %>% mcmc_plot(type = "neff_hist")

Ratios all very high.

More diagnostics

mckeon.brm3 %>% mcmc_plot(type = "combo", variable = pars)

mckeon.brm3 %>% mcmc_plot(type = "violin", variable = pars)

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

mckeon.brm3 %>% get_variables()

##  [1] "b_Intercept"           "b_SYMBIONTcrabs"       "b_SYMBIONTshrimp"     
##  [4] "b_SYMBIONTboth"        "sd_BLOCK__Intercept"   "r_BLOCK[1,Intercept]" 
##  [7] "r_BLOCK[2,Intercept]"  "r_BLOCK[3,Intercept]"  "r_BLOCK[4,Intercept]" 
## [10] "r_BLOCK[5,Intercept]"  "r_BLOCK[6,Intercept]"  "r_BLOCK[7,Intercept]" 
## [13] "r_BLOCK[8,Intercept]"  "r_BLOCK[9,Intercept]"  "r_BLOCK[10,Intercept]"
## [16] "prior_Intercept"       "prior_b"               "prior_sd_BLOCK"       
## [19] "lp__"                  "accept_stat__"         "stepsize__"           
## [22] "treedepth__"           "n_leapfrog__"          "divergent__"          
## [25] "energy__"

pars <- mckeon.brm3 %>% get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") %>% na.omit()

mckeon.brm3$fit %>%
  stan_trace(pars = pars)

The chains appear well mixed and very similar

stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

mckeon.brm3$fit %>%
  stan_ac(pars = pars)

There is no evidence of auto-correlation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

mckeon.brm3$fit %>% stan_rhat()

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

mckeon.brm3$fit %>% stan_ess()

Ratios all very high.

mckeon.brm3$fit %>%
  stan_dens(separate_chains = TRUE, pars = pars)

ggmcmc

The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

## mckeon.ggs <- mckeon.brm3 %>% ggs(burnin = FALSE, inc_warmup = FALSE)
## mckeon.ggs %>% ggs_traceplot()

The chains appear well mixed and very similar

gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

## ggs_autocorrelation(mckeon.ggs)

There is no evidence of auto-correlation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

## ggs_Rhat(mckeon.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

## ggs_effective(mckeon.ggs)

Ratios all very high.

More diagnostics

## ggs_crosscorrelation(mckeon.ggs)

## ggs_grb(mckeon.ggs)

Model validation

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_data
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

mckeon.brm3 %>% pp_check(type = "dens_overlay", nsamples = 100)

The model draws appear to be consistent with the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

mckeon.brm3 %>% pp_check(type = "error_scatter_avg")

This is not really interpretable

intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

mckeon.brm3 %>% pp_check(group = "BIRD", type = "intervals")

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(mckeon.brm2)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- mckeon.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
mckeon.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = mckeon$PREDATION,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(mckeon.resids, quantreg = TRUE)

Conclusions:

the simulated residuals do not suggest any issues with the residuals
there is no evidence of a lack of fit.

Partial effects plots

brms

conditional_effects

mckeon.brm3 %>%
  conditional_effects() %>%
  plot(points = TRUE)

ggpredict

mckeon.brm3 %>%
  ggpredict() %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

## $SYMBIONT

ggemmeans

mckeon.brm3 %>%
  ggemmeans(~SYMBIONT) %>%
  plot(add.data = TRUE, jitter = c(0.5, 0))

fitted_draws

## Partial residuals in binomial models are too confusing for the average viewer
## as they will yeild values that are not exactly 0 or 1 and this seems wrong.
## Partial.obs <- mckeon.brm3$data %>%
##     mutate(Pred = predict(mckeon.brm3, re.form=NA)[,'Estimate'],
##            Resid = resid(mckeon.brm3)[,'Estimate'],
##            Obs = Pred + Resid)

mckeon.brm3 %>%
  epred_draws(newdata = mckeon, re_formula = NA) %>%
  median_hdci() %>%
  ggplot(aes(x = SYMBIONT, y = .epred)) +
  geom_pointrange(aes(ymin = .lower, ymax = .upper)) +
  geom_line() +
  ## geom_point(data = Partial.obs,  aes(y = Obs,  x = SYMBIONT),
  ##            position = position_nudge(x = 0.1, y = 0), color = "red") +
  geom_point(
    data = mckeon, aes(y = PREDATION, x = SYMBIONT), alpha = 0.2,
    position = position_jitter(width = 0.2, height = 0)
  )

Model investigation

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

mckeon.brm3 %>% summary()

##  Family: binomial 
##   Links: mu = logit 
## Formula: PREDATION | trials(1) ~ SYMBIONT + (1 | BLOCK) 
##    Data: mckeon (Number of observations: 80) 
##   Draws: 3 chains, each with iter = 1000; warmup = 200; thin = 5;
##          total post-warmup draws = 480
## 
## Group-Level Effects: 
## ~BLOCK (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.12      1.22     1.53     6.03 1.00     2061     2049
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          2.86      1.11     0.76     5.17 1.00     2391     2259
## SYMBIONTcrabs     -1.86      0.88    -3.60    -0.16 1.00     2475     2410
## SYMBIONTshrimp    -2.31      0.90    -4.11    -0.55 1.00     2263     2305
## SYMBIONTboth      -3.16      0.97    -5.12    -1.46 1.00     2438     2328
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

the coefficients are presented on a logit scale. Whilst this is not relevant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
if we exponentiate the coefficients ($log(\frac{\rho}{1-\rho})$ -> $\frac{\rho}{1-\rho}$), they will be presented on a odds ratio scale, and thus:
- the intercept (none symbionts) will be 17.54. That is, corals without a symbiont are 17.54 times more likely to be predated on than not predated on. The odds of predation in this the absence of symbionts is 17.54:1.
- in the presence of a crab symbiont, the odds of being predated on are only 0.16 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 84%.
- in the presence of a shrimp symbiont, the odds of being predated on are only 0.1 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 90%.
- in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 96%.
if we back-transform the intercept full to the response scale (probability scale), ($log(\frac{\rho}{1-\rho})$ -> $\rho$), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1

as_draws_df (posteriors)

tidyMCMC

mckeon.brm3$fit %>%
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )

Conclusions:

see above

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

mckeon.brm3 %>% get_variables()

##  [1] "b_Intercept"           "b_SYMBIONTcrabs"       "b_SYMBIONTshrimp"     
##  [4] "b_SYMBIONTboth"        "sd_BLOCK__Intercept"   "r_BLOCK[1,Intercept]" 
##  [7] "r_BLOCK[2,Intercept]"  "r_BLOCK[3,Intercept]"  "r_BLOCK[4,Intercept]" 
## [10] "r_BLOCK[5,Intercept]"  "r_BLOCK[6,Intercept]"  "r_BLOCK[7,Intercept]" 
## [13] "r_BLOCK[8,Intercept]"  "r_BLOCK[9,Intercept]"  "r_BLOCK[10,Intercept]"
## [16] "prior_Intercept"       "prior_b"               "prior_sd_BLOCK"       
## [19] "lp__"                  "accept_stat__"         "stepsize__"           
## [22] "treedepth__"           "n_leapfrog__"          "divergent__"          
## [25] "energy__"

mckeon.draw <- mckeon.brm3 %>%
  gather_draws(`b.Intercept.*|b_SYMBIONT.*`, regex = TRUE)
mckeon.draw

We can then summarise this

mckeon.draw %>% median_hdci()

## On a odd ratio scale
mckeon.draw %>%
  mutate(.value = exp(.value)) %>%
  median_hdci()

mckeon.brm3 %>%
  gather_draws(`b_Intercept.*|b_SYMBIONT.*`, regex = TRUE) %>%
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = .variable,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)

mckeon.brm3 %>%
  gather_draws(`.Intercept.*|b_SYMBIONT.*`, regex = TRUE) %>%
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_halfeye(aes(x = .value, y = .variable)) +
  theme_classic()

bayesplot

mckeon.brm3$fit %>% plot(type = "intervals")

half-eye (ggdist)

mckeon.brm3 %>%
  gather_draws(`^b_.*`, regex = TRUE) %>%
  filter(.variable != "b_Intercept") %>%
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free")

mckeon.brm3 %>%
  gather_draws(`^b_.*`, regex = TRUE) %>%
  filter(.variable != "b_Intercept") %>%
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  geom_vline(xintercept = 0, linetype = "dashed")

density ridges (ggridges)

mckeon.brm3 %>%
  gather_draws(`^b_.*`, regex = TRUE) %>%
  filter(.variable != "b_Intercept") %>%
  ggplot() +
  geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed")

## Or in colour
mckeon.brm3 %>%
  gather_draws(`^b_.*`, regex = TRUE) %>%
  filter(.variable != "b_Intercept") %>%
  ggplot() +
  geom_density_ridges_gradient(aes(
    x = (.value),
    y = .variable,
    fill = stat(x)
  ),
  alpha = 0.4, colour = "white",
  quantile_lines = TRUE,
  quantiles = c(0.025, 0.975)
  ) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous() +
  scale_fill_viridis_c(option = "C")

## Fractional scale
mckeon.brm3 %>%
  gather_draws(`^b_.*`, regex = TRUE) %>%
  filter(.variable != "b_Intercept") %>%
  ggplot() +
  geom_density_ridges_gradient(aes(
    x = exp(.value),
    y = .variable,
    fill = stat(x)
  ),
  alpha = 0.4, colour = "white",
  quantile_lines = TRUE,
  quantiles = c(0.025, 0.975)
  ) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous(trans = scales::log2_trans()) +
  scale_fill_viridis_c(option = "C")

tidy_draws

mckeon.brm3 %>% tidy_draws()

spread_draws

mckeon.brm3 %>% spread_draws(`.*Intercept.*|b_SYMBIONT.*`, regex = TRUE)

posterior_samples

mckeon.brm3 %>%
  posterior_samples() %>%
  as_tibble()

$R^2$

mckeon.brm3 %>%
  bayes_R2(re.form = NA, summary = FALSE) %>%
  median_hdci()

mckeon.brm3 %>%
  bayes_R2(re.form = ~ (1 | BLOCK), summary = FALSE) %>%
  median_hdci()

## if we had random intercept/slope
mckeon.brm4 %>%
  bayes_R2(re.form = ~ (SYMBIONT | BLOCK), summary = FALSE) %>%
  median_hdci()

Further analyses

In addition to comparing each of the symbiont types against the control group of no symbionts, it might be interesting to investigate whether there are any differences between the predation protection provided by crabs and shrimp, as well as whether having both crabs and shrimp symbionts is different to only a single symbiont type.

These contrasts can be explored via specific contrasts.

SYMBIONT	Crab vs Shrimp	One vs Both	None vs Symbiont
none	0	0	1
crab	1	1/2	-1/3
shrimp	-1	1/2	-1/3
both	0	-1	-1/3

brms

Pairwise comparisons

mckeon.brm3 %>%
  emmeans(~SYMBIONT, type = "response") %>%
  pairs()

##  contrast       odds.ratio lower.HPD upper.HPD
##  none / crabs         6.31    0.5424     28.23
##  none / shrimp        9.98    0.2400     47.26
##  none / both         22.98    1.7479    126.33
##  crabs / shrimp       1.59    0.0673      7.56
##  crabs / both         3.59    0.2594     19.26
##  shrimp / both        2.33    0.1018     11.96
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95

## On the fractional scale
mckeon.em <- mckeon.brm3 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  pairs() %>%
  gather_emmeans_draws() %>%
  mutate(
    Eff = exp(.value),
    PEff = 100 * (Eff - 1)
  ) # ,
# Prob = plogis(.value))
mckeon.em %>% head()

mckeon.em %>%
  group_by(contrast) %>%
  dplyr::select(contrast, Eff) %>%
  median_hdi()

mckeon.em %>%
  group_by(contrast) %>%
  summarize(Prob = sum(Eff > 1) / n())

## On a probability scale
mckeon.em <- mckeon.brm3 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  regrid() %>%
  pairs() %>%
  gather_emmeans_draws() %>%
  mutate(Eff = .value) # ,
mckeon.em %>% head()

mckeon.em %>%
  group_by(contrast) %>%
  dplyr::select(contrast, Eff) %>%
  median_hdi()

## Cell means
mckeon.em <- emmeans(mckeon.brm3, ~SYMBIONT, type = "link") %>%
  gather_emmeans_draws()
mckeon.em %>%
  mutate(P = plogis(.value)) %>%
  median_hdci(P)

Specific contrasts

cmat <- cbind(
  crab_vs_shrimp = c(0, 1, -1, 0),
  one_vs_both = c(0, 1 / 2, 1 / 2, -1),
  symbiont = c(1, -1 / 3, -1 / 3, -1 / 3)
)

mckeon.em <- mckeon.brm3 %>%
  emmeans(~SYMBIONT, type = "link") %>%
  contrast(method = list(cmat)) %>%
  gather_emmeans_draws() %>%
  mutate(Fit = exp(.value))

mckeon.em %>% median_hdci(Fit)

mckeon.em %>%
  group_by(contrast) %>%
  median_hdi(Fit)

mckeon.em %>%
  group_by(contrast) %>%
  summarize(sum(Fit > 1) / n())

newdata <- emmeans(mckeon.brm3, ~SYMBIONT, type = "response") %>% as.data.frame()
head(newdata)

ggplot(newdata, aes(y = prob, x = SYMBIONT)) +
  geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD))

R^2

mckeon.brm3 %>% bayes_R2(re.form = NA)

##     Estimate Est.Error       Q2.5     Q97.5
## R2 0.1496247 0.0780289 0.01267212 0.2869891

mckeon.brm3 %>%
  bayes_R2(re.form = NA, summary = FALSE) %>%
  median_hdci()

mckeon.brm3 %>%
  bayes_R2(re.form = ~ (1 | BLOCK), summary = FALSE) %>%
  median_hdci()

## for random intercept/slope model
mckeon.brm4 %>%
  bayes_R2(re.form = ~ (SYMBIONT | BLOCK), summary = FALSE) %>%
  median_hdci()

References

Mckeon, Seabird, Adrian Stier, Shelby Mcilroy, and Benjamin Bolker. 2012. “Multiple Defender Effects: Synergistic Coral Defense by Mutualist Crustaceans.” Oecologia 169 (February): 1095–1103. https://doi.org/10.1007/s00442-012-2275-2.

Bayesian GLMM example 4

Murray Logan

04 July, 2022

Preparations

Scenario

Read in the data

Exploratory data analysis

Fit the model

rstanarm

Using default priors

Assessing priors

Defining priors

Plotting prior and posterior

brms

Using default priors

Defining priors

Plotting prior and posterior

Random intercept/slope model

Compare models

MCMC sampling diagnostics

brms

bayesplot

stan plots

ggmcmc

Model validation

brms

pp check

DHARMa residuals

Partial effects plots

brms

conditional_effects

ggpredict

ggemmeans

fitted_draws

Model investigation

brms

summary

as_draws_df (posteriors)

tidyMCMC

gather_draws

bayesplot

half-eye (ggdist)

density ridges (ggridges)

tidy_draws

spread_draws

posterior_samples

\(R^2\)

Further analyses

brms

Pairwise comparisons

Specific contrasts

R^2

References